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Two-Channel  Linear-Predictive  Spectral  Analysis; 

Program  for  the  HP  9845  Desk  Calculator 

Introduction 

Spectral  analysis  of  short  data  segments  by  the  standard  FFT  procedure  is  not  a 
viable  approach;  unstable  and/or  coarse  estimates  of  the  spectra  result.  An  at¬ 
tractive  technique  in  this  case  is  linear-predictive  spectral  analysis,  both  for  the 
single-channel  as  well  as  the  multiple-channel  cases.  See  references  1-9,  particularly 
references  7-9  which  derive  and  give  Fortran  programs  for  a  multiple-channel 
linear-predictive  spectral  analysis  technique  that  is  a  generalization  of  Burg’s 
technique  for  the  single-channel  case  (reference  1). 

The  purpose  of  this  report  is  twofold:  first,  we  translate  the  Fortran  program  in 
reference  9  into  Basic  for  use  on  the  Hewlett-Packard  HP  9845  Desk  Calculator, 
and  in  the  process,  also  make  some  minor  improvements  and  modifications  to  the 
format  and  printout  statements.  We  also  limit  consideration  to  the  two-channel  case 
and  thereby  take  advantage  of  some  simplifications  in  computing  possible  for  this 
special  case.  Second,  we  apply  the  program  to  a  pair  of  stationary  processes,  one  of 
which  has  pure  tones  that  are  not  present  in  the  other  process.  In  this  manner,  we 
point  out  a  possibly  deleterious  effect  on  the  auto-spectral  estimates  and  the 
coherence  estimate,  and  indicate  a  method  for  circumventing  some  of  the  difficulty. 
As  a  byproduct,  a  philosophy  for  multichannel  spectral  analysis  is  suggested  for 
further  consideration. 
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Use  of  Program  For  Spectral  Analysis 

In  appendix  A,  the  listing  for  the  two-channel  linear-predictive  spectral  analysis 
technique  is  presented.  Inputs  required  of  the  user  are  the  following: 

N  Number  of  data  points  in  each  process 
Pmax  Maximum  order  of  predictive  filter  to  be  considered 
Nfft  Size  of  FFT  to  be  used  in  spectral  computation.  (1) 

In  addition  to  these  integer  inputs,  the  user  must  modify  the  subroutine 
SUB  Data  (N,X(*))  to  accommodate  and  read  in  his  particular  two-channel  data. 
All  data  are  presumed  real. 

The  program  computes  the  (sample)  means  of  each  of  the  two  processes  and 
subtracts  the  means  from  the  data.  (Some  possible  ramifications  of  this  procedure 
are  considered  in  reference  6,  appendix  B;  in  addition,  the  effect  of  choosing  too 
small  an  FFT  size,  Nfft,  is  discussed  in  reference  6.)  Next,  the  covariance  matrix  (at 
zero  delay)  of  the  input  data  is  computed,  and  the  Akaike  Information  Criterion 
(AIC,  reference  8,  pages  42-44)  is  evaluated  and  used  to  select  the  integer 

Pbest  Best  order  of  predictive  filter  to  use.  (2) 

The  forward  and  backward  partial  correlation  coefficients  (references  7-9)  are 
evaluated  through  order  Pmax,  as  well  as  the  forward  predictive  filter  coefficients 
for  Pbest.  The  normalized  correlation  matrices  are  computed  through  Pmax  (ex¬ 
trapolated  values  beyond  Pbest)  and  the  spectral  density  matrix  is  computed  (via  an 
FFT)  from  zero  to  Nyquist  frequency,  fN  =  (2A)-1,  where  A  is  the  time-sampling 
increment  of  the  processes.  A  partial  check  on  the  adequacy  of  the  FFT  size,  Nfft,  is 
afforded  by  a  printout  of  the  areas  under  the  spectral  estimates  and  comparison 
with  the  (sample)  covariances  of  the  input  data.  Finally,  the  inverse  FFT  of  the 
spectral  estimate  gives  the  aliased  normalized  correlation  matrices;  the  motivation 
and  equations  for  this  approach  are  given  in  reference  9. 

A  sample  printout  for  a  short  data  sequence  (20  data  points  in  each  process)  is 
given  after  the  program  listing  in  appendix  A,  as  a  test  or  check  case  on  a  user- 
written  program.  Also,  plots  of  the  corresponding  auto-spectral  estimates  and  the 
coherence  estimates  are  given  there  for  completeness,  although  this  example  has  no 
real  physical  significance. 

Timing  Results 

Execution  times  for  the  five  major  subroutines, 

Pcc  Partial  correlation  coefficients 

Pfc  Predictive  filter  coefficients 

Peftf  Predictive  error  filter  transfer  function 

Sdm  Spectral  density  matrix 

Acm  Aliased  correlation  matrices,  (3) 
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are  given  in  tables  1-5  below,  for  the  HP  9845B  Desk  Calculator  equipped  with  the 
Fast  Processor  Upgrade  Kit.  Only  those  variables  utilized  in  each  subroutine  are 
considered  in  these  tables,  since  execution  time  is  independent  of  the  other 
variables;  for  example,  the  execution  time  of  subroutine  Pcc  does  not  depend  on 
Pbest. 


Table  1.  Execution  Times  for  Subroutine  Pcc 


N 

Pmax 

Seconds 

20 

6 

1.9 

50 

10 

5.4 

100 

5 

4.7 

100 

10 

9.5 

100 

15 

14.1 

1000 

47 

404.2 

Table  2.  Execution  Times  for  Subroutine  Pfc 


Pmax 

Pbest 

Seconds 

5 

1 

.09 

10 

1 

.15 

6 

4 

.24 

15 

5 

.62 

15 

11 

1.41 

47 

12 

4.06 

Table  3.  Execution  Times  for  Subroutine  Peftf 


Pbest 

Nfft 

Seconds 

4 

256 

17.5 

11 

256 

17.8 

1 

512 

32.0 

5 

512 

32.8 

1 

1024 

63.9 

11 

1024 

66.6 

Table  4.  Execution  Times  for  Subroutine  Sdm 


Nfft 

Seconds 

256 

8.9 

512 

17.7 

1024 

35.3 
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Table  5.  Execution  Times  for  Subroutine  Acm 


Nfft 

Seconds 

256 

9.8 

512 

18.6 

1024 

37.7 

From  these  tables,  we  are  able  to  extract  the  following  fairly  accurate  rules  of 
thumb:  the  execution  time  of 

Pcc  is  linearly  dependent  on  N  and  Pmax 
Pfc  is  linearly  dependent  on  Pmax  and  Pbest 

Peftf  is  linearly  dependent  on  Nfft,  but  is  essentially  independent  of  Pbest 
Sdm  is  linearly  dependent  on  Nfft 

Acm  is  linearly  dependent  on  Nfft.  (4) 

These  rules  allow  extrapolation  to  other  cases  of  interest  to  the  user.  The  execution 
times  of  the  FFT  itself  are  given  in  table  6. 

Table  6.  Execution  Times  for  Subroutine  FftlO 


Nfft 

Seconds 

128 

2.6 

256 

4.5 

512 

8.4 

1024 

17.1 

If  the  user  is  interested  only  in  obtaining  the  predictive  filter  coefficients  (for 
example,  to  do  time  domain  prediction  and  signal  processing),  these  results  are 
available  immediately  after  execution  of  subroutine  Pfc.  There  is  then  no  need  to 
resort  to  the  frequency  domain  routines  that  follow  Pfc;  in  this  manner,  execution 
time  and  storage  can  be  significantly  reduced.  An  additional  reduction  in  execution 
time  is  available  by  declaring  all  the  loop  counters  in  a  subroutine  to  be  INTEGER. 
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Application  to  Processes  With  Tones 

Our  first  example  is  the  two-channel  case  given  numerically  by  the  sample  values 
in  reference  7,  page  17;  reference  8,  page  K-12;  and  in  reference  9,  page  D-18.  The 
analytic  expression  for  the  autoregression  is 

x,(k)  =  .85  x, (k- 1 )  -  .75  x2(k- 1 )  +  w,(k) 

x;(k)  =  .65  x ,(k- 1 )  +  .55  x2(k- 1 )  +  w,(k)  (5) 

where  {w,(k)}  and  {w\(k)}  are  uniformly  distributed,  independent  white  noise 
processes  with  zero  means  and  variances  1/12.  General  filter  and  spectral  relations 
for  moving-average  and  autoregressive  processes  are  given  in  appendix  B;  these 
general  relations  are  then  specialized  to  this  particular  numerical  example.  It  is 
shown  in  (B-31)  et  seq.  that  the  auto  spectrum  of  process  {x,(k)}  has  four  poles  and 
three  zeros  in  the  finite  z-plane,  even  though  the  two-channel  recursion,  (5),  is  only 
first-order  regressive. 

Generally,  for  a  two-channel  P-th  order  regression  and  independent  white  ex¬ 
citations  (i.e. ,  Ek  =  0  for  k  >  P.Fk  =  and  Q(z)  =  A1  in  (B- 18)),  the  auto-  and 

cross-spectra  of  the  processes  each  possess  4P  poles  and  3P  zeros  in  the  finite  z- 
plane  (of  which  P  zeros  occur  at  the  origin).  This  is  in  contrast  to  the  single-channel 
case,  where  2P  poles  (and  a  P-th  order  zero  only  at  the  origin)  can  occur.  This  in¬ 
creased  generality  can  be  anticipated  by  the  observation  that  whereas  a  single¬ 
channel  approximation  requires  estimation  of  only  P  parameters,  an  M-channel 
approximation  requires  estimation  of  M:P  parameters  (4P  for  the  two-channel  case 
M=2).  Of  course,  for  a  fixed  number,  N,  of  data  points  from  each  process,  the 
estimation  of  an  increased  number  of  parameters  can  only  be  done  with  increased 
variance;  this  is  a  manifestation  of  the  tradeoff  between  resolution  and  stability  that 
accompanies  all  spectral  analysis  techniques. 

The  first-order  forward  partial  correlation  coefficient  for  two-channel  process  (5) 
is 

’.85  -.75* 

A* 11  (true)  =  ,  (6) 

’  L.65  .55 . 

and  all  other  higher-order  coefficients  are  zero.  The  exact  auto  spectrum  of  the  first 
process,  {xt(k)},  is  shown  (in  dB)  in  figure  1A;  the  auto  spectrum  of  the  second 
process,  {x-,(k)},  is  shown  in  figure  IB;  the  magnitude-squared  coherence  is 
displayed  in  figure  1C;  and  the  argument  of  the  complex  coherence  or  cross  spec¬ 
trum  is  depicted  in  figure  ID.  There  is  seen  to  be  a  strong  narrowband  component  at 
approximately  one-fourth  of  the  Nyquist  frequency  fN  =  (2A)-1,  where  A  is  the  time¬ 
sampling  increment  for  the  two-channel  process  (5).  This  leads  to  a  peak  magnitude- 
squared  coherence  value  of  .999013  at  2fA  =  .2459. 

The  results  of  applying  the  two-channel  spectral  analysis  program  in  appendix  A 
to  the  numerical  data  cited  above,  with  N  =  100,  are  shown  in  figure  2,  where  the 
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four  parts  of  this  figure  correspond  directly  to  those  of  figure  1 .  Pbest  turns  out  to 
be  equal  to  the  correct  value  1,  and  the  spectral  estimates  are  all  quite  good.  In  fact, 
the  estimated  magnitude-squared  coherence  reaches  a  peak  value  of  .999745  versus 
the  true  value  of  .999013. 


The  covariance  matrix  of  the  process  generated  by  (5)  is  (reference  8,  page  18, 
after  scaling  by  variance  1/12) 


R0(true)  = 


'2.095 

0.405 


0.405  ' 
1.804 


(7) 


The  corresponding  matrix  estimate  yielded  by  the  program  here,  based  on  the 
particular  N  =  100  data  values  cited  above,  is 


R 


O 


4.62 
.  .916 


.916 
3.80  . 


forN  =  100; 


(8) 


these  values  are  approximately  2.2  times  larger  than  (7),  due  to  the  fact  that  (5)  is  a 
narrowband  process  and  the  particular  100  pairs  of  samples  used  in  the  spectral 
estimates  happen  to  lie  on  a  local  peak  of  the  instantaneous  waveforms.  Although 
the  local  estimates  of  the  absolute  power  levels  are  off  considerably,  the  estimate  of 
the  forward  partial  correlation  coefficient  is  very  good;  we  find,  instead  of  (6), 


A<»  = 


.872 

-.770' 

.634 

.560. 

for  N  =  100. 


(9) 


Next  we  add  a  pure  tone  only  to  the  second  process  {x,(k)}  at  a  frequency  equal  to 
0.6  of  the  Nyquist  frequency,  i.e.,  at  0.6fN.  The  power  in  the  tone  is  1/512,  i.e.,  32.9 
dB  below  the  average  power,  3.80,  in  this  particular  segment  of  autoregressive 
process  (x2(k)};  see  (8).  The  resultant  spectral  estimates  are  shown  in  figure  3;  they 
are  virtually  identical  to  figure  2.  The  only  inadequacy  of  figure  3  is  that  the 
autospectral  estimate  in  figure  3B  gives  no  indication  of  the  added  tone;  of  course, 
there  should  ideally  be  no  indication  of  the  tone  in  figure  3A  for  the  auto  spectrum 
of  (x,(k)}.  The  value  of  Pbest  was  again  1,  as  determined  by  the  AIC. 

When  the  tonal  power  in  the  second  process  is  increased  to  -26.9  dB,  Pbest  in¬ 
creases  to  4  (see  figure  4)  and  there  are  humps  in  both  auto-spectral  estimates  near 
the  tone  frequency  0.6fN.  The  coherence  estimates  (magnitude  and  argument)  are 
significantly  perturbed  in  a  considerable  neighborhood  of  0.6fN;  this  broad 
frequency-perturbation  width  is  due  to  a  small  value  of  Pbest  having  been  selected 
by  the  AIC. 

Increasing  the  tonal  power  to  -20.8  dB  results  in  the  estimates  depicted  in  figure 
5.  Now  there  is  a  considerable  indication  of  the  tonal  power  in  figure  5B;  however, 
there  is  also  an  undesirable  indication  in  figure  5A  at  frequency  0.6fN  in  the  auto¬ 
spectral  estimate  for  process  {xt(k)}.  This  “feed-across”  is  due  to  the  fact  that  we 
are  working  with  only  N  =  100  data  samples  of  each  process;  with  this  small  a  data 
set,  the  “best”  two-channel  linear-prediction  is  misled  into  an  erroneous  indication. 
It  is  important  to  observe  at  this  point  that  any  auto-spectral  estimate  based  on 
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samples  of  process  { x,(k)}  alone  would  not  give  this  tonal  indication,  since  the  tone 
is  not  present  in  this  process. 

The  coherence  estimates  in  figure  5  fare  no  better,  even  though  Pbest  =  8  now.  A 
large  magnitude-squared  coherence  value  of  0.85  is  yielded  at  frequency  0.6fN.  The 
progression  towards  poorer  behavior  is  also  present  in  figure  6,  which  employs  a 
tonai  power  of  -14.8  dB  relative  to  the  sample  power  in  {x2(k)}.  Now  the  undesired 
peak  magnitude-squared  coherence  estimate  is  0.9.  A  tonal  power  of  -8.8  dB  (figure 
7)  yields  a  near-unity  magnitude-squared  coherence  estimate  at  0.6fN,  and  a  very 
substantial  tonal  indication  in  the  auto  spectrum  of  {x,(k)},  figure  7A. 

The  situation  is  markedly  improved  if  more  data  samples  are  available.  When  N 
is  increased  to  1000,  and  data  are  generated  via  (5)  as  before,  the  sample  covariance 
for  the  particular  data  set  generated  is 

2.60  .514' 

Ro  =  for  N  =  1000,  (10) 

L  .514  2.27  . 

for  no  tone  present.  When  a  tone  is  added  to  process  {x2(k)},  with  strength  -24.6  dB 
relative  to  the  sample  power,  2.27,  of  the  second  process,  the  resultant  spectral 
estimates  are  as  displayed  in  figure  8.  There  is  a  slight  hump  at  0.6fN  in  figure  8B, 
and  a  near-zero  coherence  estimate  at  this  frequency.  Recall  that  the  ideal 
characteristics  would  be  identical  to  figure  1  except  for  an  impulse  in  figure  8B  at 
0.6fN  and  a  very  sharp  null  in  the  magnitude-squared  coherence  at  0.6fN. 

The  results  in  figure  8  were  achieved  by  taking  Pmax  =  8,  for  which  the  AIC 
indicated  Pbest  =  8  for  this  particular  data  set.  However,  the  AIC  is  a  very  flat 
function  of  filter  order  P  in  this  range,  and  it  is  difficult  to  justify  a  particular  value 
of  P  as  “best”.  Some  additional  information  about  the  autoregressive  portion  of 
the  observed  process,  such  as  a  limit  on  P,  could  be  useful;  for  example,  when  we 
specified  Pmax  as  1,  the  results  were  very  similar  to  figure  1.  There  was  virtually  no 
indication  of  the  tone  in  any  of  the  spectral  estimates,  even  though  it  was  in  the 
{ x2(k)}  data  at  a  relative  level  of  -24.6  dB  with  respect  to  the  sample  power,  2.27,  of 
the  second  autoregressive  component.  In  fact,  the  estimated  first  partial  correlation 
coefficient  was 

‘.8543  -.7394' 

A<()  =  for  N=  1000,  (11) 

[.6578  .5415. 

which  is  very  close  to  the  true  value,  (6). 

Results  for  the  spectral  estimates  when  the  tonal  power  is  increased  to  -18.6  dB, 
-12.6  dB,  -6.6  dB,  and  -0.6  dB  are  given  in  figures  9,  10,  11,  and  12,  respectively, 
all  corresponding  to  Pmax  =  8  and  Pbest  =  8.  Even  for  the  nearly  0  dB  case  in 
figure  12,  there  is  virtually  no  indication  in  auto-spectral  estimate  12A  of  the 
strong  tonal  in  process  { x2(k) },  figure  12B.  The  magnitude-squared  coherence 
estimate  in  figure  12C  appears  to  have  developed  a  couple  of  zeros  and  poles  near 
the  frequency  f  =  0.6fN,  where  the  strong  tone  is  located;  recall  that  we  have 
4P  =  32  poles  available  in  the  approximation  for  Pbest  =  8.  Typically,  it  has  been 
observed  that  a  strong  tonal  present  in  only  one  process  manifests  itself  in  the 
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coherence  estimate  as  a  sharp  spike  at  the  tone  frequency.  The  change  in  argument 
in  figure  12D  in  the  neighborhood  of  this  frequency  can  serve  as  an  indicator  of  the 
number  of  poles  and  zeros  clustered  there. 

When  Pmax  was  increased  to  47,  the  A1C  yielded  Pbest  =  25  for  these  last  four 
figures.  However,  the  spectral  estimates  for  Pbest  =  25  proved  to  be  too  spiky  and 
erratic.  Also  the  selection  of  Pbest  at  25  is  rather  tenuous,  as  figure  13  indicates;  this 
is  a  plot  of  the  A1C  versus  filter  order  P  in  the  range  (1, 47).  Although  the  absolute 
minimum  occurs  at  P  =  25,  there  are  significant  drops  in  the  curve  at  P  =  4,  6,  and 
8.  Selection  of  P  at  one  of  these  significant  drops  appears  to  be  a  promising  ap¬ 
proach,  instead  of  using  the  absolute  minimum  of  the  curve.  In  addition,  the 
flatness  of  the  curve  is  brought  out  by  observing  that  the  range  of  values  of  AIC  is 
limited  to  (-4.80,  -4.73)  for  P  in  the  wide  range  from  10  to  47.  Thus,  the  local  minor 
drops  and  rises  in  the  AIC  curves  are  not  significant;  selection  of  values  of  P 
corresponding  to  significant  decreases  seems  to  be  a  viable  approach. 


Figure  13.  Akaike  Information  Criterion  for  N  =  1000,  Tone  Power  =  -12.6  dB 


The  last  example  we  consider  is  a  two-channel  process  of  N  =  64  data  points, 
composed  of  several  tones,  some  of  which  are  at  common  frequencies,  and  some  of 
which  are  not;  this  example  was  supplied  by  S.  L.  Marple  (reference  10).  In  par¬ 
ticular,  process  {xt(k)}  has  two  strong  tones  at  f  =  0.4fN  and  0.5 fN,  and  a  weaker 
tone  (-20  dB)  at  f  =  0.2fN,  in  addition  to  some  low  level,  colored  background  noise. 
The  other  process  has  two  strong  tones  at  f  =  0.4fN  and  0.8fN,  and  a  weaker  tone 
(-20  dB)  at  f  =  0.2fN.  Thus  the  tonal  frequencies  common  to  both  processes  are 
0.2fN  and  0.4fN,  whereas  the  uncommon  frequencies  are  0.5fN  and  0.8fN.  The  two 
auto-spectral  estimates  of  each  process  (obtained  via  the  single-channel,  forward- 
backward  averaging  technique  of  reference  4)  are  displayed  in  figure  14  for 
prediction  length  P  =  12  (24  poles  for  each  spectral  estimate).  There  is,  of  course, 
no  cross-feed  at  frequencies  0.5fN  and  0.8fN. 

The  spectral  estimates  of  the  same  two-channel  data  (obtained  via  the  program  in 
appendix  A  which  includes  coherence  estimation)  are  given  in  figure  15.  The  value 
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of  P  used  was  6,  which  allows  for  24  poles  in  each  spectral  estimate.  The  low 
number  of  data  points,  N  =  64,  now  allows  some  undesired  cross-feed  in  figures 
15A  and  15B  at  f  =  0.8fN  and  0.5fN,  respectively.  This  also  shows  up  in  the 
magnitude-squared  coherence  estimate  as  two  very  sharp  spikes  at  these  two 
frequencies,  whereas  the  true  coherence  is  zero  at  these  two  frequencies.  This 
limited  capability  of  the  multi-channel  linear  predictive  technique  can  be  improved 
by  utilizing  larger  data  sets;  N  =  64  is  too  small  a  data  size  to  accomplish  a  high 
quality  result  for  a  data  set  such  as  this  with  strong  interfering  tones. 


14A.  Auto  Spectrum  of  First  Process 


Figure  14.  Auto-Spectral  Estimates  for  Multitone  Example,  N  =  64.  P  =  12 


,t' 

t  a 
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Discussion  and  Conclusions 

A  program  for  two-channel  auto-  and  cross-spectral  estimation  has  been 
presented  and  illustrated  for  cases  including  interfering  tones.  If  the  number  of 
available  points,  N,  is  too  small,  some  misleading  estimates  may  be  obtained 
because  of  cross-feed  between  the  finite  lengths  of  data  from  each  channel.  This 
cross-feed  manifests  itself  as  narrow  spurious  spikes  in  the  spectral  and  coherence 
estimates.  Choice  of  the  appropriate  value  of  filter  order,  P,  is  possible  by  ob¬ 
serving  where  the  AIC  undergoes  significant  negative  jumps,  rather  than  by  using 
the  absolute  minimum  of  the  curve.  It  can  also  be  illuminating  to  overlay  plots  made 
with  two  (or  more)  different  values  of  P,  thereby  obtaining  different  degrees  of 
resolution  and  stability  from  the  same  data  set.  The  recursive  nature  of  the  linear- 
predictive  approach  makes  this  practice  easy  to  achieve. 

A  more  fundamental  observation  about  spectral  estimation  in  general  is  now 
developed.  Suppose  we  are  given  finite  data  records  of  three  stationary  processes 
x(t),  y(t),  and  z(t),  and  we  wish  to  estimate  all  the  auto  spectra  and  cross  spectra 
involved.  The  Blackman  and  Tukey  and  weighted-FFT  approaches  evaluate  the 
auto  spectrum  of  each  process  separately.  Thus,  the  spectrum  of  v(t)  is  estimated 
without  interference  from  y(t)  and  z(t);  the  availability  of  the  data  records  for  y(t) 
and  z(t)  plays  no  part  in  the  eventual  auto-spectral  estimate  for  x(t).  Additionally, 
the  cross-spectral  estimate  for  processes  x(t)  and  y(t)  is  independent  of  the  available 
data  on  the  z(t)  process.  Finally,  the  coherence  estimate  between  two  processes  is 
independent  of  any  additional  data  records  for  other  (statistically  related)  processes. 

On  the  other  hand,  the  generalization  (in  references  7-9)  of  Burg’s  single-channel 
linear-predictive  spectral  analysis  approach  to  the  multichannel  case  gives  auto- 
spectral  estimates  of  the  x(t)  process  that  are  dependent  on  the  available  values  of 
y(t)  and  z(t).  Also,  the  cross-spectral  estimate  between  x(t)  and  y(t)  is  dependent  on 
the  particular  z(t)  data  available.  This  procedure  can  be  poor  for  short  data  lengths 
if,  for  example,  y(t)  contains  a  strong  tone  at  fc  that  is  not  present  in  x(t)  or  z(t). 
Thus,  estimates  of  spectra  Gxx(f),  Gxy(f),  and  Gxz(f)  all  contain  tonal  indications  at 
f0  that  should  not  be  there.  These  spurious  tonal  indications  are  due  to  cross-feed 
between  the  available  finite  data  segments  of  the  various  processes. 

This  raises  the  following  questions: 

•  Should  the  estimate  of  Gxx(f)  be  determined  only  from  the  available  x(t) 
data  record  ? 

•  Should  the  estimate  of  G,v(f)  be  determined  only  from  the  available  x(t) 
and  y(t)  data  records  ? 

•  If  coherence  Cxy(f0)  =  0,  why  use  y(t)  to  estimate  Gxx(f0)  ? 

•  If  coherence  Cxy(f0)  =  1,  why  use  the  completely  statistically  dependent 
y(t)  data  to  estimate  Gxx(f0)  ? 

This  philosophy  of  discarding  “irrelevant”  data  would  be  consistent  with  the 
Blackman  and  Tukey  and  FFT  approaches.  Carrying  this  philosophy  on,  we  are  led 
to  the  following:  estimate  Gxx(f)  solely  from  the  x(t)  data  by  some  good  single- 
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channel  linear-predictive  technique,  such  as  forward-and-backward  averaging, 
coupled  with  an  efficient  way  of  inverting  the  relevant  matrices  (e.g.,  references  4 
and  5).  Then  estimate  cross  spectrum  Gxy(f)  or  coherence  C(f)  directly,  by  some  (yet 
unknown)  linear  predictive  technique  whose  sole  goal  is  linear  prediction  of  x(t) 
from  y(t)  and  vice  versa,  with  no  interest  in  or  diversion  from  simultaneous 
estimation  of  Gxx(f)  or  Gyy(f).  By  this  means,  we  can  concentrate  on  extracting  all 
the  relevant  cross-spectral  information  with  maximum  stability  and  resolution. 
Other  cross  spectra  of  interest  between  particular  pairs  of  available  processes  can  be 
similarly  obtained,  one  at  a  time.  This  procedure  is  currently  under  investigation. 


ii 
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Appendix  A 

Program  For  Two-Channel  Linear-Predictive  Spectral  Analysis 

The  program  listing  below  in  Basic  for  the  HP  9845B  Desk  Calculator  is  a 
translation  and  update  of  that  given  in  references  7-9.  A  complete  breakdown  and 
explanation  of  the  components  and  subroutines  of  the  program  are  given  in 
reference  8,  and  in  reference  9,  appendix  D. 

Inputs  required  of  the  user  are  the  integers  listed  in  lines  20,  30,  40;  they  are 

N  Number  of  data  points  in  each  process; 

Pmax  Maximum  order  of  predictive  filter  to  be  considered; 

Nfft  Size  of  FFT  to  be  used  in  spectral  analysis. 

In  addition,  the  user  must  modify  subroutine  SUB  Data(N,X(*))  in  lines  5430-5490 
at  the  end  of  the  program  to  read  in  his  own  particular  two-channel  data  sets.  Pbest 
can  be  forced  to  equal  Pmax  by  setting  Fac  =  0  in  SUB  Pcc. 

An  explanation  of  the  program  output  is  given  under  equation  (1)  of  the  main  text 
of  this  report.  A  sample  printout  for  a  short  (N  =  20)  data  sequence  that  can  be 
used  as  a  check  case  on  the  program  is  presented  following  the  listing  below.  Sample 
plots  of  the  auto-spectral  estimates  and  the  coherence  estimates  conclude  the  ap¬ 
pendix. 
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10  1  rwO-C  HANNEL  LINEAR  PREDICTIVE  I  F  EC  TF  AL  “flHL.fi:.,  rP  “ 5  0 1  ,  * 
20  N  =  20  •  NUMBER  OF  BATA  POINT:  If)  EACH  PROCESS 

30  F'fii  a  =  6  i  MAXIMUM  ORDER  OF  PREDICTIVE  FILTER 

40  Nf  f  r  =250  !  SIZE  OF  FFT 

50  OPTION  BASE  1 

60  R E D I M  Y '  t ,  II  1  ,  Z 1  2  ,  H  ,  m p,  ■  Pm  a  ,2,2  > ,  E p  Fi  j  .  2 , 2 

70  REDIM  Rn  •  Pm  a  ,  2, 2  1 ,  H>  'V  0:  F'iuj  i « , X  1 1  ■  Nt  f  t  ■ ,  Y  1  1  ■  Nf  t  •  .  "  1  2  ■  (if  *' » 

SO  REDIM  Y  1  2 1  Nt' f  t  i ,  X21  •  Hfft  '• .  Y21  •  Nt  ft',  X22'  Nf  ft  ■ ,  .  22  N*  f  »  1 

90  DIM  Y  <  2 , 1 0OO  > , Z  C  2 ,  1 000  > , A p  t  25 . 2 , 2  ■ . B P ' 2 5 , 2 . 2  ■ 

1  00  DIM  Rn  •:  25 , 2 ,  2  ' ,  A 1  c  C  0 :  25  '  ,  X 1  1  <  1 024  > ,  V  1  1  ■  1 024  ■ . X 1 2 ■  1 024  • 

110  DIM  Y 12  •  1  024  >  ,  X2 1  <:  1024  > ,  Y2 1  •  1024  f ,  X22  ■  1024',  ,22'  1024' 

120  DIM  Ave ■ 2 > , Ubs  s  t ' 2, 2 > , U < 2 , 2  :■ , V <  2 , 2  > , U t • 2 , 2  ■ , v t • 2 . 2  • , A ■ 2 , 2  ■ 
130  DIM  B<  2, 2’> .  P'  2, 2)  ,  Wa<  2, 2>  ,  Wb-  2.2).  Nc  '2.2',  Md-  2.2*.  Wee  2, 2  • 

140  PRINT  "NUMBER  OF  DATA  POINTS  IN  EACH  PROCESS  M  =";N 

150  PRINT  "MAXIMUM  ORDER  OF  PREDICTIVE  FILTER  Pma  =";Pma 
160  PRINT  "SIZE  OF  FFT  Nt'ft  =  "  ;  Nf  ft 

170  PRINT 

ISO  CALL  Dit i(H,Y,t)) 

190  PRINT  "PROCESS  NUMBER  1" 

200  FOR  1=1  TO  N 

210  PRINT  YU,  I  >5 

220  NEXT  I 

290  PRINT  L I M <  1  > 

240  FRINT  "PROCESS  NUMBER  2" 

250  FOP  1=1  TO  N 

260  PRINT  Y  <  2 ,  I  > ; 

270  NEXT  I 

2S0  PRINT  l I N < 2  > 

290  CALL  Pc  c  ( N ,  Pm  a  <: ,  V'  *  • ,  Z<  *  ) .  A“«  >  *  > ,  Ha'  *  ,  nt  •*•,  Uc  Md  ,  W« 

<  *  ) ,  A  i  c  >  *  ) ,  P  b  e  1 1  ,  U  b  *  s  t  <  *  > ,  U  t  <  *  > ,  v  i  •  *  ’ ,  A  ■  *  ■■ ,  B  ■  *  ■ ,  A  p  •  *  .Ip"  *>  • 

300  PRINT  "MEANS  OF  INPUT  DATA  <.Ave->:" 

310  PRINT  Rued) 

320  PRINT  A',e<2)  ' 

330  PRINT  L I N( 1 > 

340  PRINT  "COVARIANCE  MATRIX  OF  INPUT  DATA  ' P  )  1 " , R 
350  PRINT  " AKAIKE  INFORMATION  CRITERION:" 

360  PRINT  “  P  AicfPi" 

370  IMAGE  3  D ,  4  <  4  X ,  M .  9DE ) 

3S0  FOP  P  =  0  TO  Pm*,; 

390  PP I  NT  USING  370; P, At c ( Pj 

4O0  NEXT  P 

410  PRINT  LINel) 

420  PRINT  "Pbest  *“ ; Pbest 

430  PRINT  LINO) 

440  PRINT  "Ubeit : “ , Ubes> <  *  ) 

450  PRINT  "FOPWAPD  PARTIAL  CORRELATION  COEFFICIENTS:  ' 

460  PRINT  "  P  Apt,  1,1)  ApC  2,  1  '  Apt, 

Ap  f  2 , 2 ) " 


470 

FOR  P»1  TO  Pma< 

4  S  0 

PRINT  USING  370JP, ApCP,  1, 1), Ap' P,2, 1 >, Ap- P, I 

, , 2 ■ . Ap  P , 2 

490 

NEXT  P 

500 

PRINT  LINO) 
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510  PRINT  "BROWARD  PARTIAL  CORRELATION  COEFFICIENT:" 

520  PRINT  "  P  Bp'  1,1;  Bp' 2, 1 <  Bp  1 , 2 > 

Bp '  2 , 2  *  " 

500  FOR  P=1  TO  Pm*x 

540  PR  I  NT  US  I  NG  370  ;  P ,  Bp '  P  ,  1  ,  1  >  ,  Bp  • P ,  2 ,  1  1  ,  Bp 1  P  ,  1  ,  2  ■  ,  Bp  ■  P  ,  2 ,  2  > 

550  NEXT  P 

SCO  PRINT  LIN''1> 

570  IF  Pfcc i t  =0  THEN  890 

58 0  CALL  Pf  c  (. Proa* ,  Pbe si ,  RC  *  > ,  Ac > :  *  ■ .  Bp<:  *■  • ,  u«i  *  • .  Wb-  *  • ,  Wc  ■  *  ■ ,  »lo •  *  ■ .  Pru  *  ■ ,  R 1  1 ,  Pi! 

,R12> 

5*0  PRINT  "FORWARD  PREDICTIVE  FILTER  COEFFICIENT'S  FOR  Pb«*t:“ 

600  PRINT  "  p  Ap  <1,1)  Hp '0,1'  Ap<l,2> 

Ap<2, 2 > " 

610  FOR  P«1  TO  Pb«  st 

620  PRINT  USING  370 ;  P ,  Ap P,  1 ,  1  >  ,  Ap  i P ,  2.  1  > ,  Ap •  P,  1 , 2  ;  ,  Ac ■  P  ,  2  .  2  ; 

630  NEXT  P 

640  PRINT  L I N < l > 

650  PRINT  "NORMALIZED  CORRELATION  MATRICES  cRr. ; :  • 

660  PRINT  “DELAY  AUTO  11  CF0SS21  CROSS  12 

AUT022" 

670  PRINT  USING  370;  0,  R  ■’  I ,  1  >  ,  R  <  2 ,  1  ,  P  <  1,2  *  ,  fi '.  2 , 2  • 

680  FOR  P=1  TO  Prnax 

6*0  PRINT  USING  3  70 ;  P  ,  Rn  <  P ,  1  ,  1  > ,  Rn  <  P , 2 ,  1  ■ ,  Rr,  •  P ,  1,2  • ,  Pn  •  P ,  2 , 2  • 

700  NEXT  P 

710  PRINT  LIN( 1 ) 

720  CALL  P#f  t  f  (Pbest ,  Nf  f  t ,  Ap<*  • ,  XI 1  c  *  » ,  Y  1 1  <  *  • ,  XI  2'  *  •  .  V 1  2  ■  ♦  1  ,  X2 1  •.  *  ' ,  Y2 1  »  *  > ,  X22  ' 
* ) , Y22< ♦  >  > 

730  CALL  Sdm  <  Nf  f  »  ,Ub«*»  <1  *  >  ,  U*<  *  > ,  Wb<  *  ,  Uc  ■  * '• ,  W<3<  •*>  ,  U«  •  *  '.111'  *  > ,  V  1  1  '.  *  .> ,  X  1  2  <  *  1  , 
Y 1 2 '  *  ■ ,  X2 1  >  * ) ,  Y2 1  i  *  > ,  X22 ,  Y22  •  *  > ,  5 1 1 ,  $22 ,  S 1 2  ) 

740  PRINT  "SPECTRAL  DENSITY  MATRIX  AND  COHERENCE,  FROM  ZERO  FREQUENCY  'BIN  1C 

750  PRINT  "  BIN  AUTOU  AUT022  RE'.  CROSS  12  •>  IN'CR0SS12>  MAG  SC  CGH 

ARGUMENT" 

760  IMAGE  3D,5<M.6DE, 1X),M.6DE 
770  FOR  1=1  TO  30 

780  L-I 

790  IF  I< 16  THEN  840 

300  IF  I >16  THEN  830 

813  PRINT  "***" 

820  GOTO  850 

830  L-I+Nf f t '2-29 

840  PRINT  USING  760;L,X1 1 CL>,X22' L> ,X1£' L  •  , Y12' L> .  .1  1  ■ L  ■ .  . 22<L> 

850  NEXT  I 

860  PRINT  LINC1) 

870  PRINT  "TRAPEZOIDAL  SUMS  OF  SPECTRA:" 

880  PRINT  S11,S22,S12 

890  PRINT 

900  PRINT  “COVARIANCES  OF  INPUT  DATA:" 

910  PRINT  Rll.R22.R12 

920  PRINT  L IN< 1 > 

930  CALL  AcmCNfft ,  XI  1  <*>  ,  X12<  *'» ,  Y12C  *>  ,  X21  <  *  > .  Y21  ■  -  '  ,  '22'  -  • .  XI  1ml ,  X22ml ,  XI  ImO, 

X22m0> 

940  N 1 “Nf  f t ♦ 1 

950  N2*Nf  f t  s 2 

960  N22-N2+2 
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970  PRINT  "ML  IH3ED  NORMALIZED  C  OPPELH  T  I  ON  MATRICES:  " 

980  PRINT  "  DEL  Mi'  MUTOll  CROSS Ol 

MU  TO  22 “ 

990  PRINT  USING  370; O, X 1 1 <  N02  1 , XO 1 '  1  ',..21-  1  ■  ,  XOO ■  ’■Li  ■ 

1  Ciiao  FOR  1  =  1  TO  07 
1010  L  =  I 

1000  IF  I  16  THEN  1070 
1030  IF  116  THEN  1060 
1040  PRINT  “*■**“ 

1O50  GOTO  1030 
1O60  L-I+N2-29 

1070  PRINT  USING  3”0 ;  L  , : :  1  1  >  NOO  +  L  '  ,  XO  1  '  N  1  -  i»  ■ , :  0  1  ■  1  '  00  II 

1000  NEXT  I 

1030  PRINT  USING  370; NO- 1 , X 1 lml , KOI  ■  NOO > , KOI ■ NO  ■ . XOOn,  1 
llOO  PRINT  USING  370;  NO,  XI  lm0,  XOl  (.110+  1  ' ,  XOl  •  NO"- 1  • .  ..OOdiO 
1110  PRINT  LINO' 

1100  PRINTER  IS  0 

1130  PRINT  “ HUT 0  SPECTRAL  DENSITIES  IN  DEI  ' 

1140  PLOTTER  IS  “GRAPHICS" 

1150  GRMPHICS 

1160  SCALE  0, NO, -5,0 

1170  GRID  N2/4, 1 

1180  PENUP 

1190  FOR  1=0  TO  N2 

1000  '"LOT  I  ,  LGT<  XI  1<  I  +  l  >  ) 

1010  NEXT  I 

1200  PENUP 

1230  FOR  1=0  TO  N2 

1240  PLOT  I, LGT<J<22<I41  >) 

1250  NEXT  I 
1260  PENUP 
1270  DUMP  GRAPHICS 
1280  PRINT  LIN<3> 

1290  PRINT  "MAGNITUDE  SQUARED  COHERENCE  HnD  ARGUMENT" 

1300  PLOTTER  IS  “GRAPHICS” 

1310  SCALE  0.N2.0, 1 
1320  GRID  N2'4, .25 
1330  PENUP 
1340  FOR  1*0  TO  N2 
1350  PLOT  I,YU<I  +  1> 

1360  NEXT  I 
1370  PENUP 

1380  SCALE  O , N2 , -P I , P I 
1390  FOR  1=0  TO  N2 
1400  PLOT  I» Y 22(1+1) 

1410  NEXT  I 
1420  PENUP 
1430  DUMP  GRAPHICS 
1440  PRINT  L I N  <  4 ) 

1450  PRINTER  IS  16 
1460  END 
1470  I 
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1430  :UB  Pc  c  '  N ,  Phi  a  ,  Y  •  *  '  ,  2  (  *  1  ,  A"<  ■  •  •  .  W  *  ■  ,  Wb  ■  »  1  .  W  *  ■  »  .  L< 

*  ,  A i  ;  ■'  *  ' ,  Pbeat ,  Ubsat  ■  *  1 ,  Ui  •.  *  > ,  V l  ■  *  • ,  A-  * ' ,  B-  *  • ,  Ap  *  .  Bp •  * 

1440  Ia*INT( 1 . 5*SQP(N>  ' 

1500  IF  Pm  ax .  =  I  a  THEN  1520 

1510  PRINT  "Pm a. .  *"  ;  Pma.<:  "  IS  TOO  LARGE  FOP  N  =  " ;  N  ;  -  SEAR 
1520  I a=M I N< I  a,  Pma*.  > 

1530  Fac=8-N  '  Fac=0  WOULD  FORCE  Fb«  =  *  EOohl  10  Rr.,a 
1540  MAT  Aoe=RSUM'Y> 

1550  NAT  Ave»Aoe- (Ml 
1560  A 1 =Ave  <  1  ) 

1570  A2  =  Aoe  <  2 > 

15S0  FOR  I«1  TO  N 
1590  Y<1,  I>»YU,  I)-A1 

16O0  Y<2,  n»Y(2,  I  >-A2 

1610  NEXT  I 
1620  MAT  Z*Y 

1630  CALL  Aut o < 2 , N- 1 , Y ( * > , Wc t * > > 

1640  CALL  Aut ot 1 , 1 , Y< * > , Wd( * > i 
1650  CALL  Aut  o( N, N, Y<  *  > , Uc *  > 1 
1660  MAT  Wa=Wc  +  We 
1670  MAT  Wb=Wc  +Wd 
1680  MAT  R=*Wb  +  W» 

1690  MAT  R*R  '  ( N ) 

1700  MAT  U*R 
1710  MAT  V*R 

1720  CALL  Cr  osa  <  2 ,  N ,  Y  c  *  )  ,  2  <  *  1 ,  Wc  •  •*  >  > 

1730  At c (01-LOG( DET(U) > 

1740  A i e m i n«A i c < 0 > 

1750  Pb«st»0 
1760  MAT  Ub*st=U 
1770  FOP  P«1  TO  la 
1730  MAT  V i *  I NV (  V  ) 

1790  MAT  Ud*Vi*Wb 
1300  MAT  Wb-Wd 
1810  MAT  U 1  * I NV  <  U ) 

1820  MAT  Wd«Wa 
1830  MAT  Wa»Wd*Ui 
1840  MAT  Wc  *Wc  * ( 2 ) 

1850  CALL  So  1  oe(  Wac  * ) ,  Wb  C  *  >  ,  Me  ( * ) ,  Wd<  *  > ,  We  ‘ '*  '■  1 
I860  MAT  A«Wc*Vi 
1870  MAT  Wd«TRH<Wc> 

1880  MAT  B«Wd»Ui 
1890  Ap(P, 1 , 1 )«A( 1 ,  1) 

1900  Ap<  P , 1 , 2  >  *A< 1 , 2  > 

1910  Ap(P,2, 1 ) »A( 2 , 1> 

1920  Ap <  P ,  2 , 2 )“A«. 2 , 2  > 

1930  8p(P, 1, 1>»B<1,  1) 

1940  Bp<P, 1,2)-B<1,2) 

1950  Bp<P,2, 1)«B<2,  1) 

1960  Bp<P| 2, 2>»B<  2,2) 

1970  MAT  W*»A*Wd 
1980  MAT  U*U-W# 

1990  MAT  W*»B*Uc 
2000  MAT  V»V-W« 
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2  0 1  0  Hi  c  >.  P  >  =  L  0  G 1  BET  >  U  1  >+F  ac  +  P 
X  0  X  0  IF  H  i  •:  ■  F'  1  =H  i  c  in  i  n  THEN  X  O  6  0 
2030  Hicnnn  =  Hic<P) 

2040  F'b«St=P 
2050  MAT  Ub«st»U 
2060  IF  P=  I  a  THEN  XXXu 
2070  L=P  ♦  1 

2080  FOR  K=N  TO  L  STEP  -1 
20*0  R 1 »  r-  •  t ,  K  • 

2100  H2*Y  <2,10 
2110  Bl*2a,K-n 
2120  B2*2 ■ 2 , K-  1  ' 

2130  Z< 1 , K  >  *B 1 -B ■  1  ,  1  <  *  A 1  -  B  <  1 , 2  ’  ■*  A  2 
2140  Z(2,lO»B2-B>2,  1  > *A 1 -B <  2 , 2  >  +  A2 

2150  VU  ,  1  ,  1  *»B1  -A<  1 , 2  >+B2 


2160  Y<2,IO«A2-A<2,  1>*B1-A<.2,2>  +  B2 

2170  NEXT  K 

2180  CFILL  Hut  oCP+2,  N,  Y<  *  > ,  Ua<  *  >  '• 

2190  CRLL  Hut o < P  + 1 , N- 1 , 2 v * > , Wb > *  •  ■ 

2200  CALL  Cross'.  P  +  2,  N,YC»), 2  <*>, Wet*') 

2210  NEXT  P 

2220  A 1  = .  5*' Ubest  >  1 , 2)+Ub#st  <  2.  1  -■  < 

2230  Ub*st < 1 , 2)*Ubest <2, 1 >*H1 
2240  3UBENB 
2250  ! 

2260  SUB  Pf  c  ■- Pma  . ,  Pbest ,  R <  *  > ,  Hp 1  +  • ,  Bp ■  +  1 ,  U  a  -  .  i’-t"  *  .  im  ■  -  ■ ,  WO',  *  1 ,  fin  v  *  ■ .  R 1  I  ,  P22 . 
R  1  2  > 

22  70  Rr,<l  ,  1  ,  t  '*Ap<l , 1 , 1 >*R<  1 ,  1  >+Hp>.  1 ,  1 ,2  -*P-  2,  1  ■ 

2280  Rr><  1 ,  1,2  1  =  flp  <  1 ,  1,  1  » *R  <  1 , 2  '>  +Ap  1 ,  1,2  ■  ~P  <  2,2  • 

2290  Rm  1,2,  l'*«Hp<.  1,2,  1>*R<1,  1  >+Ap*  1 , 2, 2  >+R  •  2.  1  • 

2300  Rn<.l,2,2'=Hp''l,2,l  >*R<1,2)+Ap'.  1,2, 2  >  ♦  P 2,2  ■ 

2310  FOR  P«2  TO  Pbest 

2320  Uc  <  1 ,  1 ) =  Ap < P , 1 , 1 >*R<  1 , 1 >+Hp' P, 1 , 2i*P' X.  i  ' 

23  30  Uc  <  1 ,  2->=Ap<  P,  1 , 1 >*P<  1 , 2>  +Ap<P,  1 , 2)  *R‘  X  ,  2  • 

2340  Uc  (2,  l  >=Ap<P,  2,  1 >*R< 1 , 1 )+Ap<P,2,2v*R. 2. 1 > 

2350  Uc <2, 2>=flp<P, 2, l>*R(l,2)+flptP,2,2)*R<2.2) 

2360  FOR  l_»l  TO  P-1 

2370  Ib»P-L 

2380  U*a,  1)»HP'.P,  l,  l;*Ep(Ib,  l.  1  >+Hp<P,  1,2  -Bp  lb.  2,  1  ■ 

2390  Ua< 1,2'  *Hp'P, 1, O  +  Bptlb, 1 , 2)+Ap<P, 1,2'  -  Bp'  lb, 2,2  ■ 

24  00  Ui(  2 ,  1  ’  *Ap  <  P  ,  2  ,  1  '■  *Bp  <  I  b ,  1  ,  1 1  +Hp*  P ,  2 , 2  ■  +  Ep  >  I  c  ,  2  .  1 

2410  Ua<2,2i*Ap<P,2,  1  *Bp <  I b,  l,2>+Ap'P,2,2  '-Bp'  Ib.2,2 

2420  U*( 1 , 1 ;*Hp<L, 1 , 1 >-Ua< 1 , l > 

2430  UtM,  2  >>flp<L,  1,2>-U*<1,2> 

2440  Ua>  2, 1  )=flp<L,  2,  1 -  Ua'  2 ,  1 ) 

2450  Uat2, 2>»flptL,  2, 2>-W*<2, 2) 
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2490  Mb'  1 , 1  '  =Bp  i P,  1  ,  1  >  *hp  >  L  i  1.1  *  ♦Bp  ■  P  ,  1 , 2  L.2.  1  • 

24  70  Wb  •  1 , 2  ■  =  Bp  ■;  P ,  1 ,  1  > ♦  Ap <  L ,  1 ,  2  >  +  Bp ■  P ,  1 , 2  •  -ftp •  >_,2,2  • 

2480  Mb'  2,  1  -■  *Bp< P ,  2,  1  ■' ♦  h p 1  L ,  1 ,  1  ■  ♦Bp'  P  ,  2 . 2  *  »«,•  •  L  >  2  .  1  • 

2490  Wb '  2,  2'>=Bp‘  P,2,  1  l,2j+Bp'P.2,2  1  »ftp  ■  L .  2  •  2 

2500  Bp  *  I  b,  1 ,  1  >=Bp< Ib.  1 .  1  >-Wb<  1 ,  1  •> 

2510  Bpc  lb,  1 , 2  >-hf><  lb,  1 ,  2'>-WbCl ,  2> 

2520  Bp<  Ib,2,  n=Bp<  Ib,2,  1  >-Wb<  2,  t  > 

2530  Bp<Ib,2,2>*Bp<.Ib,2,2>-Ub*.2,2,« 

2540  Ap<L,  1 ,  l  >«U*<.1,  1> 

2550  flpCL,  l,2>*W*a,2.' 

2560  flp(L,2,  1  )«Wa<2,  1  ,■ 

2570  flpa,2,2'>*W*i2,2> 

2580  Ud^  1 ,  1  >»Wa<,  1 ,  1  '  *Rn<  lb ,  1 ,  1  >+Hi'  1 , 2  >*Pr,.  Ib,  2.  1  ■ 

2590  Wd< 1,2>*Wai 1 , 1  >*Rn<  Ib,  1 , 2  '+Ui'  1 , 2  '♦Pr,.  Ib.  2, 2' 

2800  Ud*'2,  1  >*W*C2,  l>*Rr.Ub,  1 ,  1  )+Wv  2,2  ■  » R f. •  Ib,2,  1  • 

2810  Wd 1  2 , 2  )»U»«  2 ,  1  )*Rn(  Ib,  1 , 2  ■♦Wal  2, 2’-*Piv  Ib,  2, 2  • 

2820  MAT  Wc*Uc+Wd 

2630  NEXT  L 

2640  Rn<P, 1, l'«Wc< 1 , 1 > 

2650  RruP,  1 , 2  J«Wc  1  1 , 2  ' 

2660  Rn <  P, 2 , 1 ) *Wc <2,17 
2670  Rn<  P, 2, 2>*Wc (2,2; 

2680  NEXT  P 

2690  FOR  P*Pb«si  +  l  TO  Pm»x 
2700  MAT  Ua-ZER 
2710  FOR  L* 1  TO  Pb«st 
2720  Jb»P-L 

2730  Mb (  1 ,  1  '<  =fip ( L  i  1  ,  n*Rru  Ib,  1  ,  l  >+Ap«.L,  1 , 2  •♦Riv  Is,  2,  1  1 

2740  Mb<  1 , 2  >  *Ap '  L  ,  1  ,  1  >*Rrulb,  1 . 2  '  +Ap<  L ,  1,2  '*P.rv  Ib,  2. 2  ■ 

2750  Ub>2,  1  >=Ap'.  L,  2,  1  >*Rn<  Ib ,  1  ,  1  <  +  Ap<.  U ,  2 , 2  :•  *Prv  I  b  ,  2  .  1  • 

2760  Ub<2,2>«fip<L,2, 1 >*Pn< lb, 1 , 2 ) *ftp< L , 2 , 2  •♦Rn>  Ib, 2. 2  • 

2770  MAT  W*»U*+Ub 

2780  NEXT  L 

2790  Rn<P,l,l>«U»<l,l> 

2800  Rn<P, 1,2>»U«<1,2> 

2810  RnCP,2, 1 >-W*<2, 1 > 

2820  Rn<.P,2,2>*W*<2,2? 

2830  NEXT  P 
2840  R11*R<1,1> 

2850  R22*R<  2, 2  > 

2860  R12»R<1,2’> 

2870  T-SQRCRl 1*R22 J 

2880  R< 1 , l >»R<2, 2)»1 

2890  R<1,2)»R<2, 1>«R12'/T 

2900  FOR  P»l  TO  Pm*x 

2910  RruP,  1 ,  1  )“Rn<P,  1 ,  1 ;  -  R  1 1 

2920  Rn  <  P , 1 , 2 ) "Rn  <  P , 1 , 2 ) ^T 

2930  Rn<P, 2, 1 >»Rn( P, 2, 1 )^T 

2940  Rn<P,2,2)»Rn<P,2,2)^R22 

2930  NEXT  P 

2960  SUBEND 

2970 
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2  9  SO 

2990 

3000 

3010 

3010 

3O30 

3040 

3050 

30<?0 

3070 

30:30 

3090 

3100 

3110 

3120 
1 2  v  * ; 
3130 
3140 
3150 

3  1  60 
3170 
3130 
3190 
3200 
3210 
3220 
3230 
3140 
3250 
3260 
3270 
3230 
3290 
3300 
3310 
3320 
3330 
3340 
3350 
3360 
3370 
3380 
3390 
3400 
3410 
3420 
3430 
3440 
3450 
3460 
3470 
3480 
3490 
3500 
3510 


SUE  Pef  t  *  ■  Ft-*  *  *  ,  Mt  *  »  ,  ftp  ■  *  • ,  1  1  •  *  ■  ,  'i  1  1  •  »  ,  ;2  -  ■  .  .  i  1  ■  *  , :  :2  1  ■  *  ■  ,  i'2  1  ■  *  • , :  .22  ■  * 

'.*)■> 

>:  1 1 '  1  >  =  x  2  2  ■  1-«*1 

FOF  L= 1  TO  Pb«i> 

XI  l  v  L  ♦  l  ■  =  -Ap  i  L  .  1.1' 

X 1  2  v  L  +  1  >  =  - Ap  >  L .  1,2  1 
”21 '  L+l  ’*-Ap'.L  .  2,  1  > 

X22  •  L  +  1  >  *-Ap(  L ,  2, 2  > 

NEXT  L 

CALL  Ff  t  10vNt  f  t  ,  XI  1  v*  )  .  Y1  1  (  .  >  > 

CALL  Fit  1  O  1  N  f  f  t ,  X 1 2  ‘  *  1  ,  Y  1  2  ■  *  > 

CALL  FT  t  1  0  v  N  f  f  t  ,  X  2  1  *  .  *  .* ,  )  2  1  1  *  ■  ' 

CALL  F  f  t  1  O  v  N  f  f  t  ,  X  2  2  v  *  )  ,  Y  2  2  ( •*  >  '■ 

SUBEHD 


SUB  SdKi'  Nf  t' t  ,  Ub«*t  •  *  ' ,  Wa‘  *  • ,  Mb'.  *  .■ .  He  ■  *•  ■ .  lid’  *  .  We  •  *  ) , : :  1  1  •  ♦  > ,  V 1 1  ■:  *  > ,  X 1  2 >  •*  > ,  V 
,  X2  1  1  *  1  ,  Y2 1  1  *  1  ,  X22  t  *  1 ,  Y22 1  *>,311,  S22  ,'312' 

T  =  2  Nt't't 
SI  1 =S22*  3 1 2  =  0 
>Nfft<  2+1 
FOR  L=1  TO  J 
Wav  1 ,  1  >  =X22<  L  ■> 

Wav  1 ,  2)  =  -X12vL> 

Wa(  2 ,  1  j  =  - X 2  1  CL) 

Wa<  2 , 2) *X1 1 <L) 

Wb  v  1 ,  1  )  *  Y  2  2  C  L  > 

Wb< 1,2>»-Y12<L> 

Wb<  2 , 1 )*-Y2 1 <L ) 

Wb<2, 2>*Yl 1  CL) 

Ta*DET  (Wa)-DET  ( Wb ) 

Tb=Wa(  1 ,  1  >*WbC2, 2>+Wa<2, 2->*Wt"  1 ,  1 )  -W*'  1 . 2  •  *Wt"  2 .  1  •  -Hp  2,  P  *Wb  v  1 , 2 ) 

Ta*T^ vTa*Ta+Tb*Tb) 

MAT  We  =  TRN<  Wa) 

MAT  Wd=Ubest*Wc 
MAT  Wc=Wb*Ud 
Tb-Wc  < 1 , 2)-We  C2,  1  > 

MAT  Wc*Wa*Wd 
MAT  Wd*TRN<  Wb ) 

MAT  W«*Ubest  *Wd 
MAT  Wd*Wb*We 
MAT  Wc*Wc+Wd 

Y1 1(L)*(  Wc  < 1 , 2 1 " 2  +  T b * T b ) x > Wc  v  1  .  1 >*Uc  '2.2'' 

Y22  < L ) *FNAr  g <  We < 1 , 2) , Tb) 

XI 1 < L)*Ta*Wc v 1 , 1 > 

X22vL)=Ta*Wc(2,2) 

X12(L)*Ta*Wc (1,2) 

Y12CL)*Ta*Tb 
SI 1*S1 1+X1 1CL) 

S22*S22+X22(L) 

S12=S12+X12(L) 

NEXT  L 

SI  1*31 1-. 5*CX1 1 > 1 >+XllC J)  ) 

S22*S22-.5*(X22(  1 '  +X22  <  J ) ) 

S12»S12-. 5*<X12(  1  >+X12< J)  ) 

SUBEND 
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3520  SUB  flc hi •  m*  *' »  , 1 1 \  *  • ,  x i z,:  *  > ,  >  l *  ■  *  • ,  :■ -:2  1  ■  »  ■ ,  2 1  •  •»  • , :  221  *  • ,  :  1 1 r« l .  '22f»  l , . .  1 1  n> 0 , 

22t»0 > 

35  30  H2*Ufft  2 
3546  H21=N2+1 

3556  N28=N2+2 

3560  FOR  L»1  TO  N2 
3570  X21 >U)». 5+X1  l  <  L  > 

3580  V21  <C>».  5*M22<L  t 

3530  X21>  N2*L?».?*>:U'N22-L> 

3666  Y2 1  •  N2  +  l',*.5*X22<N22-L> 

3610  NEXT  L 

3620  CALL  Fft 10kNfTt  , X21 • *>,Y2l 1  *  -  > 

3636  Ta=l'X21 1 1 '■ 

3640  Tb*l^Y21tl> 

3650  T*Si3R<  T  a*Tb> 

3660  X1  1CN22>*:<22>  N22 '  =  1 
3670  FOP  L  =  2  TO  N2-1 
3680  XU^N2l+L>=X2UL  .*Ta 
3690  X 2 2 ' M  2 1  +  L  >  =  Y  2 1  <  L  >  ♦  T  b 
3700  NEXT  L 
3710  Xllml=X21 <N2>*T» 

3720  X22ml«Y2KN2>*Tb 
3730  Xllm0*X21<N21>*T» 

3740  X22m0=Y21 <  N2 1 J*Tb 

3750  X2ia>«.5*Xl2a>*T 

3760  V210  7»-.5*Y12UT4T 

3770  FOR  L  =  2  TO  N2 
3780  X21  iL>«X12a'«T 

3790  Y2 1 iL) ■- Y 1 2  <  L  >  *T 
3800  X2 1 C N2+L >  ■  Y2  1 <  N2+L  >»0 
3610  NEXT  L 

3820  X21<N2l  ■  =. 5*X 1 2 < N2 1 > *T 

3830.  Y2UN21)  — ,5*Y12<N21>*T 

3840  CALL  Ffl  10<NFft  ,  X2I<*>,  Y2t  (*,■> 

3850  SUEEND 
3860  ! 
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3870  SUB  Croi ; ■ Ml , Mi , H'  ♦  •  , S'  ♦  ■  .  C ■  *  ■  1 
3880  S 1 1=3 1 2  =  S2 1 =322=0 
3890  FOR  k=Ml  TO  Hi 
3900  fil=fi<l,k> 

3910  fi2=fl'2,K> 

3920  B  1  =  B  ■  1 , k -  1  > 

3930  B2  =  B ■  2 ,  k  -  1  '■ 

3940  Si  1*31 1 *fil »Bl 
3950  S  1  2  =  3 1 2+fi 1 *B2 
3930  S21=321+fii*Bl 

3970  S22=322+fii*B2 

3980  NEXT  k 
3990  Ca,l)  =  Sll 
4000  C(l,2v«S12 

4010  C  <  2 , 1 >  * 32 1 
4020  C '  2 , 2  v  =  S  2  2 
4030  SUBEND 
4O40  i 

4050  SUB  RutO'Nl,H2|ft'*^|B<*^> 

4080  311=312=522=0 

4O70  FOR  K=N 1  TO  N2 
4080  fil=fi(l,k> 

4090  fi2=fiv2,K> 

4100  S 1  1  *S 1 1 *fil *fi  l 

4110  S 1 2=S 1 2+fl 1 *fii 

4120  $22=S22+fi2*fi2 

4130  NEXT  K 

4140  B< 1,1V =311 

4150  B  <.  1 , 2  v  *  B  ( 2 ,  tv*  3 12 

4160  B  <  2 , 2  >  =522 

4170  SUBEND 

4130  i 

4190  SUB  SoU'fCW'  »  '  ,  Mb’  *  V  ,  Wc  '  *  > .  Wd'  * 
4200  T**W»<  1  ,  1  >♦«*•  2,2>*MbU  ,  1  V4Mb<  2, 

4210  Tb=DETvW«;-DET< WbV 
4220  MAT  Ud=Wc»Wb 
4230  Ue<  1 ,  1  >=U*-.2,2.' 

4240  Ue< 1 , iV'-WaC 1 , 2 > 

4250  U«<2,  1>»-U*<.2,  1  v 

4260  W«v 2, 2/*W*f 1, 1 V 
4270  MfiT  UU=W**Uc 
42S0  MfiT  Ud=Wd+U» 

4290  MfiT  Ub=Wb*<Ta; 

4300  Wb<  1 ,  l)=Wb<  1,  1  >*Tb 
4310  Ub<  2, 2^=Wb<2, 2v+Tb 
4320  MfiT  Ue=INV<UbV 
4330  MfiT  Wc=Wd*U* 

4340  SUBEND 
4350  1 


,  Ws  *  * 
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4 :<5 o  sub  Fft  io>  m,  >  1  n  *  2  io  =  i o 1 4 ,  !!*:•  integer 

4  370  DIM  CU:257> 

4  380  DATA  1 ,  .  999981  1  7528  3,  .  939*2470  1  5  39,  .  *  *  *  8  JO?  :  ;  7  *6  .  .  9*9  6  *  8  S  1  8  6  St,  .9995294175 

01  ,  .  939322  334503,  .  999077727753,  .  9  93795456205  ,  .  9984 75580573 ,  .  998  1181  12*00 

4  3 4  0  D  0  T  ft  .  9  9  7  7  2  3  0  6  6  6  4  4  ,,  9  *9  7  2  9  0  4  5  6  6  7  9 9  *  6  8  2  0  2  9  9  2  9  1  ,  .  9  9  6  3  1  2  6  1  2  1 8  ?  ,  .  9  9  5  7  6  7  4  l  4  4  6  :: 

,  .  9 9 5  i  S 4 7 2 6 6 7 2 ,  .  9 9 4 5 6 4 5 7 O 7 3 4 , . 9 9  3 9 8 6 9 7 O 0 0 2 . . 9 9  2 2  1  1  9 4  9 2  3 5 .  .  992473534599 

4  4  00  Dft Tft  ,991709753669, .990902833425, . 9900582  1  0252 .  . 98  9 1 765 09965 , . *88 2375677 3 1 

,  .987301410153,  .936308097245,  . 935277642  389 , . 9642 1 00 92  I  37 ,  . 98  3 1 0548743 1 

4410  Dft Tft  . 98 19638691 10, . 930785280403. . 979569’65685 . . 973317370720, . 977028142653 

, . 973702130039, . 974339332786, . 972939952206, . 97 l 50 3690996 , . 970031253 1 95 

4420  Dft Tft  . 963522094274, . 9669764 7 1 045 , . 965 3 94 4 4 1 698 , . 963776065795, .962121404269 

, . 960430519416, . 958703474896, . 956940335732. . 955 141 168306, . 953306840354 

4430  DftTft  . 951435028969, . 949528180593, . 947365591018. . 945607325 38 1 , . 943593458162 

,  .  94  154  4  065183, . 9 3 94  59223602 , .9  3733901  191  3,. 9  3518  35099 39. . 932992798835 

4440  DftTft  . 938766961079, . 928506088473, . 9262 102421 36, . *2387953251 1,. 92151 4039 34 2 

, .919113651690, . 9 l 667905992 1 , . 914209755784 , . 311706032005, . 309167983091 

44  58  DftTft  . 906595704515, . 90398929  3123, . 901  34884  7046. . 8  986744  65694, . 89596624  9756 

,  .  893224  301  196, . 8904  4372 3245 , . 837639620483 , . 8847*^0*8431 . . 331921264  34  3 

4  460  DftTft  . 379012226429, . 876070094195,  . 3730*497  34  13, . 370036991  109, . 86784  624  5516 

, . 863972856122, . 86086693S638 , . 3577236 1 0008 , .354557*38365, .351355193105 

4470  DftTft  .948120  344803,,  .34435  3565250,  .341554*774  37.  .338224  705555,  .  834362874986 

, . 831 469612303, . 828045845258, . 824589382785. . 321 1025149*1 , .81758481315 2 

4430  DftTft  .  814036329786, .810457198253, . 88634755 3544 , . 80 320753 1 48 1 , .79953726910S 

, . 795336904609, . 792106577308, . 738346427627, . 734556597156, . 730737223572 

4490  DftTft  . 776883465673, . 7738 1 045 336 3 , .769103337646. . 765 1 67265622 , . 76 1 202385434 

,. 757  208346506, . 75  3186799044, . 749 1 36 394523 , . 745057735441, . 74  0951  125355 

4500  DftTft  .  736816563877,  .  73265427 1 672,  . 728464 390443 . . "2424708295 1 , . 720002507961 

, . 715730825284, .711432195745, . 707 10678 1 1 37 , . 782754744457 , . 698 376249409 

4510  DftTft  .6939714608*8, .689548544737, . 68503 366777 3 , . 630600997795 , . 676092703575 

, .671558954847, .666999922  304, . 6624 1 5777590 , . 657306693297, . 653  1  72842954 

4520  DftTft  .648514401822, .643831542890, .639124444364, .634393284164, .629638238915 

, .624359438142, .620057211763, .615231590581, . 6 1 0332386276, . 6055 1 1 04 1 484 

4530  DftTft  .600616479384, .595699304492, .59075970135*. .58579785745 6, .580813958896 

, .575308191418, .570780745387, .565731810784, . 56066 1 576 1 97 , .555570233020 

4540  DftTft  .550457972937, .545324983422, . 540 1 7 1 472730 . . 534997619887, .529303624636 

, . 524589682678, .519355990166, . 5 1 4 1 027441 93, . 503830142543, . 503538383726 

4550  DATA  .498227666973, .492898192230, . 487550160148, .4:32133772079, .476799230063 

,  .  471  396736826, . 465976495768, . 460538710953, . 45503  3587126, . 44  961  1329655 

4560  DATA  .444122144570, . 4336 1 6238539 , .433033313853, .427555093430, .422000270880 

, . 416429560098, .410343171058, .405241314085, . 399624 1 99846 , . 393992040061 

4570  DATA  .338345046699, . 382683432365 , . 3770074 1 02 16 , . 371317193952, .365612997305 

, . 359395036535, . 354163525420, . 3484 1 86802 4 9 , . 342660717312. . 336889853392 

4530  DftTft  . 331 106305760, . 325310292162, . 319502030816, . 31 3681740399, . 307849640042 

, . 30238594  9319,  . 296  1  5088S244 ,  . 290284 677254 , . 28440753721  1  ,  . 2785  1  9689 385 

4590  DftTft  . 2726213554  50, .266712757475,. £607  941  17915,  . 254865659605 , . 24  8927605  74  6 

, . 242980179903, . 237023605994, . 23 1 058 1 0623 1 . . 22508391 1 360 . . 2 1 9 1 0 1 240 1 57 

4600  DftTft  . 21  31  10319916, . 2071 1  1  376192, .201  104  6  34  842. .  1 *5090322016, .  1S9063664150 

, . 183039887955, . 177004228412, . 1 70961 888768, . 1 6491 3 1 28498, . 158858143334 

4610  DftTft  . 152797135258, . 146730474455, . 1 48658239 3 33 , . 1 34588783507, . 128493110794 

, . 122410675199, . 1 16318630912, . 110222207294, . 1 04 1 2 1 6 3 387 2 , . *80 1 7 1 403296E- 1 

4620  DftTft  . 91 903956497  IE- 1 , . 8579731 23444E- 1 , .  796324 3797 1 4E- 1 , . 73564  5635997E- 1 , . 

674439 1 95637E- 1 , . 6 1 320736 3022E- 1 , . 55  19524  4 3497E- 1  .  . 490676743274E-1 

4630  DftTft  .  429332569349E-1 , . 3630722  294 1 4E- 1 , .  3067430 3 1 766E- 1 ,  . 2454 1 2285229E- 1  ,  . 

18406729 9058E-1, . 1 227 1 5382857E- 1 , . 6 1 35884649 1 5E-2 , O 
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4t'40  PEhB  Civ 
4650  M  =  1  0  2  4  -  N 
4660  FOP  I=0  TO  H  4 
4670  C  I  ♦  1  >  =C  M*  1  ♦  1  > 

4600  NEXT  I 
4600  N1*N  4 
4700  N2=N1+1 

4710  N  3  =  N  2  +  1 
4720  N4*N 1  +N3 

47  30  Log2ri=INT  <  1 . 4427*L0C<  N  l  +  .  5  • 

4740  FOP  11  =  1  TO  Log2r. 

4750  I2*2A<Loa2i"i-I  1  ) 

4760  13=2*12 

4770  I4=M'I3 

4780  FOR  15=1  TO  12 
4700  16='  15-1  '*14+1 

4800  IF  1 6  <  =  N  2  THEN  4840 
4810  N6=-C<N4-I6; 

4820  N7  =  -C  <■  I  6-N 1  > 

4830  GOTO  4860 
4840  N6  =  C  < 1 6 ) 

4850  N7  =  -C  ''  N3-  1 6  ) 

4860  FOP  17=0  TO  N-I3  STEP  13 
4870  18=17+15 

4880  18=18+12 

4800  N8  =  X  < I  8 ) -X  < 10) 

4000  H0  =  Y<  I8;-Y<10  ■> 

4010  X< I8)=X< I8;+X< 10) 

4020  Y(I8)«Yll8)+Y(!0) 

4030  X< I0)=N€*N8-N7*N0 
4040  Y( 10  )  =  N6*N0  +  N7*N8 

4050  NEXT  17 

4060  NEXT  15 

4070  NEXT  II 

4080  Il=Log2n+l 

4000  FOR  12=1  TO  10  i  2  10=1024 

5000  C  < 1 2 )  =  1 

5010  IF  1 2 >Log2n  THEN  5030 
5020  C< 1 2 ) =2  '  <  1 1 - 1 2 > 

5030  NEXT  12 
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S040 

K*  1 

5050 

FOR  I 1-1  TO 

C'  10  ' 

5060 

FOR  12-11  TO 

CO' 

STEP  C<.10< 

50  ?0 

FOR  13*12  TO 

C  <  8  ■■ 

STEP  C>  9  > 

5030 

FOR  14-13  TO 

C '  7  '■ 

STEF  CO' 

5030 

FOR  15*14  TO 

C  *.  6  > 

STEP  C<7  - 

5100 

FOR  16-15  TO 

C  1 5  ' 

STEP  CO' 

5110 

FOR  17-16  TO 

CO> 

STEP  C  >  5  • 

5120 

FOR  18-17  TO 

CO' 

STEP  C>  4  ■ 

5130 

FOR  19-18  TO 

C  O'* 

STEP  CO 

5140 

FOR  110=19  TO  CU>  STEF  CO' 

5150 

J-l  10 

5160 

IF  K  J  THEM 

5230 

5170 

fi*X\h  > 

5160 

X<K ;*X(  J  > 

5190 

X<  J)=fl 

5200 

fl-V(K) 

5210 

Y  <  K'  >  -V  <  .*  ■> 

5220 

Y< J'«P 

5230 

k-K  + 1 

5240 

NEXT  110 

5250 

NEXT  19 

5260 

NEXT  18 

5270 

.NEXT  17 

5280 

NEXT  16 

5290 

NEXT  15 

5300 

NEXT  14 

5310 

NEXT  13 

5320 

NEXT  12 

5330 

NEXT  11 

5340 

SUBEND 

5350 

! 

5360 

DEF  FNPrglX, 

Y  > 

•  PR  I  NCI  PPL  PR 

5370 

IF  X-0  THEN 

P» . 5*P1 *SGN< Y  ' 

5380 

IF  X< >0  THEN 

P-PTNs Y^X) 

5390 

IF  X<0  THEN 

P=fi+PI*< 1-2«<Y;0) ) 

5400 

RETURN  P 

5410 

FNENB 

5420 

i 

5430 

SUB  Bit  *<.  N ,  X 

<  *>  ) 

5440 

OPTION  BPSE 

1 

5450 

REDIM  X  <2 , N ) 

5460 

DPT  6  1,2, 6, 3 

»  1  1  * 

«*  »  1  »  ^  ,  5  «  ^  ^  1 

5470 

DPT P  2, 1,0, 1 

,  5 , 3  > 

0  t  3  j  5 . t , 

5480 

REPD  X<*> 

5490 

SUBEND 
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NUMBE»  OF  DATA  POINTS  IN  EACH  PROCESS  N  =  2‘) 

MAXIMUM  ORDER  OF  PREDICT  I/E  FILTER  Pi„  .  6 
SIZE  OF  FFT  NT  ft  =  256 

PFOCE53  NUMBER  1 

12  6  3  1  121453215  6  12-4 

PROCE 3S  NUMBER  2 

2  1015301262242  3  5  6  0 


MEANS  OF  INPUT  DATA  <A»«.‘: 

3.  35 
3.2 

COVARIANCE  MATRIX  OF  INPUT  DATA  < R  - : 
5.7275  .03 

.63  6.16 


AKAIKE  INFORMATION  CRITERION: 
F  A  i  c  <  P  > 

0  . 354 36 3690E +0 1 

1  . 324147843E+01 

2  . 31320954SE+01 

3  . 342425365E+01 

4  . 235777553E+01 

5  . 2485S9498E+0 1 

6  . 266 1171 85E+0 1 


Pb«st  »  4 


Ub*st : 

.  545 33344383  .226987553  32  3 

.226387553923  4. 00240808163 


FORWARD 

PARTIAL  CORRELATION  COEFFICIENTS: 

P 

fip<  1,  i) 

Apt  2 , 1 ) 

Apt  1 , 2  > 

Ap  1 2) 2  > 

1 

.  32684786 I E+00 

-. 296055028E-0 1 

.43":  5  0  £  9  5  8  E  ♦  0  0 

.  492296 1 0 lE  +  00 

2 

-.  497079412E  +  00 

. 851055338E-O1 

.  177974 150E  +  00 

- .  2750 1751 2E  +  O0 

3 

-. 1 72083 1 46E+00 

. 23733535SE+00 

.  586 121 047E-0 1 

•  15480372 3 E  +  00 

4 

-. 1 4363 1 6SSE+00 

-. 1 10686655E+OO 

.  7 3 474362 9 E+00 

-  .  1 24  7 1 5363E  +  O0 

5 

-.516644665E-01 

-.440540782E+00 

. 1 93329039E+00 

.  333560 194E  +  00 

6 

. 901306062E-01 

-. 73661 6304E-0 1 

-. 25653O155E+00 

.  734752737E  +  O0 

BACKWARD  PARTIAL  CORRELATION 

COEFFICIENTS: 

P 

Bp< 1,1) 

Bp  1 2 ,  1 ) 

Bp- 1.2) 

Bp  ( 2 ,  tL  > 

1 

.  39169448 1E  +  O0 

45 1 665794E+0O 

-. 1 397 17314E-01 

. 427449481E+00 

2 

-.657343877E+00 

178422353E+00 

. 27 14022 10E-01 

-.  19214351 2E  +  00 

3 

-.  218767381E  +  00 

685358 1 77E-8 1 

. 1 77465 355E+00 

.  129759954E  +  00 

4 

-.  128639908E  +  0O 

1 006 1 7477E+0 1 

-. 3655 J4759E-81 

-.94414 1 856E-0 1 

5 

- .  257 1 37623E-0 1 

1 94 1 24206E+O0 

-.21 9204 3O3E  +  0O 

. 1 17920510E+0O 

6 

. 182731463E+00 

28446316 3 E+00 

.  *  *;  a.  &.  y .  i  E  -  u  i 

. 161491 6  79E  +  00 
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forward  predictive  filter 

p  ftp  <  1 , 1  .> 

1  .  235726249E+00 

2  -.6O04S3360E+O0 

3  -.40436O932E+00 

4  - . 1 4363 1 688E+00 


COEFF  I C  I  Eh  TS  F OR  Pti 
Ap<2,  1  ■ 

.1741 18149E+U0 
- . 1 64936772E +00 
. 342089409E+O0 
-.  1  1 O6S6655E+00 


ftp  ■  1  .  2  ■ 

2  -  F.  5  7  0 1  5  E  +  0  0 
4  1  2  6 9  6  7  9  E  +  0  0 
2  ;  ?  7  4  F  4  2  E  ♦  0  0 
7  3  4  ~  4  7  F  2 9  E  +  0  0 


ftp' 2,2 ' 

.'5701  70572E  *00 
4424  1  33  30E*O0 
. 25 144*57  3E  +  00 
l247153i3E»00 


NORMAL 

IZED  CORRELATION 

MATRICES  <Prv>: 

DELAY 

AUTO  1  1 

C  R  0  S  S  2 1 

CROSS  12 

0 

. 100000000E+01 

.  1  39734996E  *-00 

■  1  3 97  349 9 6  E  +  00 

1 

. 339669762E+0O 

. 402437203E-01 

.  49525094  6  E ♦ 0  0 

2 

-. 234751305E+00 

. 91248091OE-01 

.511 09C660E  +-00 

3 

-. 372945922E+00 

. 1 95203 1 45E+08 

.  11710711 3E+0O 

4 

223941701E-01 

. 87623866 >£"01 

.  1 907927 16E+0Q 

5 

. 396697 3 60E+00 

89679 1093E "01 

.  1  1270 1 30  3E+00 

6 

.  3  3  0  3  S  S  4  O  0  E  +  O  0 

836262522E-01 

-.•=>51S6912’7E-01 

RUT022 

. 1 00000000E+0 1 
. 438307 047E+O0 
. 640872368E-O 1 
. 663036683E-O1 
. 1 00976458E+08 
. 122691975E+00 
. 1 8833678 1 E-0 1 


SPECTRAL  DENSITY  MATRIX  AMD  COHERENCE,  FROM  ZERO  FREQUENCY  'BIN  1>: 


BIN 

AUTO  1  1 

AUT022 

RE  v  CROSS 1 2  > 

IN* CROSS  12 » 

MhG  so  c  oh 

ARGUMENT 

1 

.  563234E-Q1 

. 1 304S0E+00 

.  848704E-01 

»  0  0  0  0  0  OE  ♦  0 1 

. 9801 1 7  E + 0  0 

. 0O00U0E+01 

2 

. 562462E-01 

. 1 30295E+O0 

. S471 1 3E-01 

- .  ;:54S07E-02 

*  ? S 0 0 C-4E  +  U0 

-. 30O463E-01 

3 

. 560 152E-0 1 

. 1 29740E+0O 

. 342  359E-0 1 

-.50€?5?E-0I 

. 979904E+O0 

-. 601 10SE-O1 

4 

. 356327E-0 1 

. 1 23822E+00 

. 334493E-81 

- . 754340E-O2 

.  9  7  9  6  3  4  E  ♦  0  O 

-. 902094E-O1 

5 

. 55 1 025E-0 1 

.  1 27549E  +  00 

.  32 3603E-0 1 

- . 996 1 O5E-02 

. 97  9250E+00 

-.  120360E+00 

6 

.54430 1  E-0 1 

. 1 25934E+00 

. 30981  IE-01 

-.  12287  3E-01 

. 97S746E+O0 

-.  1 5053 1 E  +  00 

7 

. 5  36226E-0 1 

.  1  2  3994E  +  0O 

. 793277E-01 

-.  1 45083E-0 1 

.9781 13E+00 

-.  130892E+0O 

3 

. 526890E-0 1 

.  121751 E+00 

. 7741 95E-0 I 

166075E-01 

.  977  342E+00 

-.21 131 1E+0O 

9 

. 5 1 6400E-0 1 

. 1 1 9230E+00 

. 752793E-Q1 

-. 185704E-01 

.  97  64 1 9E  +  00 

-.241 85SE+00 

10 

. 504882E-O 1 

. 1 16460E+00 

.  729330E-8 1 

- .  203S57E-0 1 

. 97533 1 E+00 

- , 272557E+08 

1  1 

.  492476E-0 1 

. 113474E+00 

.  704OS3E-0 1 

22044SE-01 

. 974O59E+00 

-. 303429E+00 

12 

. 4793  35E-0 1 

. 1 10309E+00 

.  677369E-0 1 

- . 2  35426E-0 1 

. 972584E+0O 

- . 334499E+00 

13 

. 465623E-01 

. 1 07002E+00 

.  649485E-01 

-. 248774E-0! 

. 970333E+00 

-. 365795E+00 

14 

. 451507E-01 

. 1 83593E+00 

. 62075 1 E-01 

- . 260503E-0 1 

. 968928E+00 

-.397' 45E+00 

15 

. 437 1 53E-0 1 

.  1001 2 1E  +  00 

.  59 1 473E-0 1 

-. 270676E-O1 

. 966668E+00 

-. 4291S1E+0O 

♦  *  # 

116 

. 495654E-02 

.527756E-02 

. 180088E-82 

. 366375E-02 

. 638529E+00 

. 1  1  1448E  +  01 

117 

.4641 26E-02 

. 492382E-02 

. 1 90232E-02 
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Appendix  B 

General  Filter  and  Spectral  Relations 

For  each  integer  k,  let  Hk  be  a  rectangular  matrix  of  (complex)  constants,  to  be 
called  the  filter  impulse  response  at  delay  kA,  where  A  is  the  time  sampling  in¬ 
crement.  The  number  of  filter  outputs  need  not  equal  the  number  of  filter  inputs. 
For  multichannel  filter  excitation  Wn  at  time  nA,  the  filter  output  at  time  nA  is  given 
by  the  discrete  convolution 


-  E 

k 


Hk  Wn-k 


(B-l) 


where  the  summation  extends  over  all  nonzero  summands. 

For  a  stationary  excitation,  let  the  correlation  matrix  of  complex  input  process 
{ Wn}  at  delay  lb  be 


W  W9  =  P  (non-diagonal  matrix) , 

*  (B-2) 

where  the  overbar  denotes  an  ensemble  average,  and  the  superscript  H  denotes  a 
conjugate  transpose.  The  z-transform  of  input  correlation  sequence  {P^  }  is 


a(z)  =  A  £2'*  P  , 
l  * 

and  the  spectrum  of  input  process  { Wn}  is,  for  real  frequency  f, 

Q(f)  =  Q.(exp(i2irfA))  =  A  exp(-i2iTfAJt)  P. 

I 

=  A  £  exp(-iu£)  P  , 

S,  * 

where  we  let 


(B-3) 


(B-4) 


u  =  2tt£A 


(B-5) 


The  correlation  mahix  of  the  filter  output  process  { Xn } ,  at  delay  lb,  is  given  by 
using  (B-l)  and  (B-2): 


R'Z  = 


X  X‘ 


n-j, 


k,m 


H.  W  , 
k  n-k 


w”  ,  h”  sj]  h.  p 

n-y.-m  m  k 

k,m 


hh 

2,+m-k  in 


(B-6) 


The  z-transform  of  output  correlation  sequence  { }is,  by  use  of  (B-6), 
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G(z)  =  A  y  z'V  =4r  z'1  r  H.  P„  ,  ^ 

T  1  l  kTm  k  ^+m-k  m 

•  <3-7, 


Now  define  the  z-transform  of  filter  sequence  {Hk}  by 

H(z)  =  £  z'k  H, 
k  K 

and  define  the  quantity 

H  H 

H  (z)  =  [H ( z) ]  ,  including  conjugation  of  z. 


Then 


«"&) -["&)]  "■[?  -k  »k]H  -  E  ^ 


(B-8) 


(B-9) 


(B-10) 


We  now  employ  (B-8),  (B-3),  and  (B-10)  in  (B-7)  to  obtain 
G( z)  =  H(z)  d(z)  HH  (p-) 

The  spectrum  of  output  process  {Xn}  is  then,  at  real  frequency  f, 

G(f)  =  G(exp(i2irfA) )  =  A  £  exp(-i2TTfA&)  R* 


(B- 11) 


=  H(f)  Q(f)  Af)  . 


where  we  employed  (B-ll),  (B-4),  and  set 


H(f)  =  H(exp(i2TrfA))=  £  exp(-i27rfAk)  ft. 

k  K 


(B-12) 


(B-13) 


We  also  employed  the  property  that 


HH(f)  =  [H(f)]H  =  £  exp(i2rrfAk)  H?  , 


(B-14) 
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or 


HH(exp(i27tfA))  =  H^f) 


(B- 15) 


Finally,  from  convolution  (B-l),  we  obtain  the  z-transform  of  output  data 
sequence  {Xn}  as 


X(z)  e  £  z‘n  Xn  =  H(z)  W(z)  , 

n  (B-16) 

where  we  used  (B-8)  and  defined 

W(z)  =  £  z-n  Wn  . 

n  (B- 17) 

The  major  results  thus  far  are  given  by  (B-l),  (B-16),  (B- 11),  and  (B-12)  for  a 
general  filter  and  excitation. 

ARM  A  Process 

For  a  multichannel,  autoregressive,  moving-average  process,  the  recursion  is 
given  by 


X  =  F  E,  X  .  +  £  F,  W  , 
n  k  n-k  “  k  n-k 

k  k 


The  z-transform  of  this  equation  is 


where 


X(z)  =  E(z)  X(z)  +  F(z)  W(z)  , 
E(z)  i  L  z‘k 


Ek  ,  F(z) 


z'k  F, 


Then  we  can  solve  (B- 1 9)  for  X  (z)  as 

X(z)  =  [I  -  E(z)]"1  F(z)  W(z)  =  H(z)  W(z) 
where  the  transfer  function  from  input  to  output  is 

H(z)  E  (I  -  E(z)]"1  F(z) 


(B-18) 


(B-l  9) 


(B-20) 


(B-21) 


(B-22) 


in  terms  of  the  parameters  of  recursion  (B-18).  But  now  (B-22)  is  in  the  framework 
of  the  presentation  above;  namely,  the  spectrum  of  output  process  {Xn}  is,  from 
(B-12), 

G(f)  =  H(f)  Q(f)  ^(f)  ,  (B-23) 
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H(f)  =  [I  -  E(f)  ]  1  F(f)  , 

E(f)  =  E(exp(i2trfA))  =  £  exp(-i2trfAk)  E 

k  K 

F(f)  =  F(exp(i2irfA) )  =  2Z  exp(-i2irfAk)  F, 

k  k  (B-24) 


Example 

As  an  example  of  (B-18),  consider  the  multichannel  first-order  recursion 

X  =  E,  X  .  +  W 
n  1  n-1  n 

with  the  input  spectrum  for  { Wn}, 

Q(f)  =  A  I  for  all  f  . 

This  is  a  white  process,  uncorrelated  from  channel  to  channel.  Then 


(B-25) 


(B-26) 


E(z)  =  z"1  ,  F(z)  =  I  , 

H(z)  =  (I  -  z_1  E^"1  ,  (B-27) 

G(f)  =  A  H(f)  H^f)  • 


Specialization  to  Two-Channel  Process 

We  further  specialize  example  (B-25)  to  the  two-input,  two-output  channel 
process  characterized  by  coefficient  matrix 


E 


1 


b 

d 


(complex  coefficients) . 


Then  (B-27)  and  (B-22)  yield  transfer  function 


(B-28) 
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H(z) 


1  -  z"1  (a  +  d)  +  z~2(ad  -  be) 


1  -  z-1  d 


z'1  c 


z-1  b 


1  -  z  *  a 


(B-29) 


There  follows 


HH(1/z*)  =  [H(1/z*)]H 


1  -  z(a  +  d)*  +  z  (ad  -  be)’ 


1  -  zd*  zc* 


zb* 


1  -  za* 


(B-30) 


and 


fi(z)  -  4  N(z)  «**(!/!*)  5 


guU) 

*21(l) 


«12U) 

*22W 


1  *  |b|2  *  I d[ 2  -  z'1  d  -  zd* 
-ab*  -  cd*  +  z"1  c  +  zb* 


-a*  b  -  c*d  +  z-1  b  +  zc* 

1  +  [ a [ ^  +  |c|2  -  z~l  a  -  za* 


(B-31) 

where 

D  «  [l  -  z_1(a  +  d)  +  z'2(ad  -  be)]  [l  -  z(a  +  d)*  *  z2(ad  -  be)*] 

(B-32) 


Inspection  of  denominator  D  in  (B-32)  reveals  that  G(z)  has  poles  (i.e.,  all  the 
elements  of  matrix  G  (z)  have  poles)  at 
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z  = 


a  +  d  +.  \  (a  +  d)  -  4(ad  -  be) 
2 


(B-33) 


and  at 


1 


z 


★ 

00 


exp(i  arg(zM)) 


(B-34) 


That  is,  even  though  recursion  (B-25)  is  limited  to  first-order  regressive  and  white 
independent  excitations  (see  (B-26)),  G(z)  has  four  poles  in  the  finite  z-plane. 
Generally,  if  H(z)  has  a  pole  at  z^,  then  W(l/z*)  has  a  pole  at  1/z* ;  so  if 


z^  =  r  exp(i0) ,  then  l/z£  =  exp(iQ)/r 


(B-35) 


has  the  same  argument. 

In  addition,  element  g,,(z)  in  (B-31)  has  a  zero  at  00  and  three  zeros  in  the  finite 
z-plane,  at 

I  .  0  and  i  ■  1  *  N2  ♦  ld|2  VCl  *  M2  -  |d|2)2-4|d|2 
0  0  2d* 

(B-36) 

The  product  of  the  latter  two  zero  locations  is  d/d*  =  exp(i2arg(d)).  Thus,  the  auto 
spectrum  of  process  1  has  three  zeros  and  four  poles,  even  though  the  multichannel 
recursion,  (B-25),  is  first-order  regressive.  Similar  behavior  is  true  of  the  auto 
spectrum  of  the  second  process,  as  well  as  the  cross  spectrum  between  the  two 
processes. 

The  magnitude-squared  coherence  for  this  example  is,  from  (B-3 1 )  and  (B-5), 
g12(exp(iu))  g21(exp(iu)) 

g11(exp(iu))  g22(exp(iu)) 

i  2 

1-  a*b  -  c*d  +  b  exp(-iu)  +  c*  exp(iu)] _ 

[|l  -  d  exp(-iu) | 2  +  | b| 2J  1  -  a  exp(-  iu) | 2  +  | c ] 2 J 

(B-37) 


This  has  four  zeros  and  four  poles  in  the  finite  z-plane. 
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Numerical  Example 

We  now  specialize  (B-28)  to  example  (6)  in  the  main  text: 


- 

— 

- 

a  b 

.85 

-.75 

c  d 

.65 

.55 

. 

(B-38) 


The  four  poles  of  the  spectrum  G( z)  are  located  at 


zm  =  ,7+i  .6819  ,  1/z*  =  .7330  +  i  .7140 


(B-39) 


and  the  zeros  of  gn(z)  are 

Zq  =  0  ,  3.0646,  .3263,  00 


(B-40) 


The  zeros  of  g12(z)  are  located  instead  at 

zQ  =  0  ,  .8802,  -1.3109,  « 


(B-41) 


The  magnitude-squared  coherence  simplifies  to 

1.0634  -  .056  cos (u)  -  .975  cosf2u) 
[1.865  -  1.1  cos(u)]  [2.145  -  1.7  cos(u)] 


(B-42) 


This  example  was  used  frequently  in  the  main  body  of  this  report.  The  peak  value  is 
.9990128  at  u  =  .772564,  or  2fA  =  .245915. 
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